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Abstract 

In a number of previous papers [20l |23l |2TJ [H |T2j |T3l [TTl |T0], local (coarse grid) multi¬ 
scale model reduction techniques are developed using a Generalized Multiscale Finite Element 
Method. In these approaches, multiscale basis functions are constructed using local snapshot 
spaces, where a snapshot space is a large space that represents the solution behavior in a coarse 
block. In a number of applications (e.g., those discussed in the paper), one may have a sparsity 
in the snapshot space for an appropriate choice of a snapshot space. More precisely, the solution 
may only involve a portion of the snapshot space. In this case, one can use sparsity techniques 
(16101a EH HU El]) to identify multiscale basis functions. In this paper, we consider two such 
sparse local multiscale model reduction approaches. 

In the first approach (which is used for parameter-dependent multiscale PDEs), we use local 
minimization techniques, such as sparse POD, to identify multiscale basis functions, which are 
sparse in the snapshot space. These minimization techniques use li minimization to find local 
multiscale basis functions, which are further used for finding the solution. In the second approach 
(which is used for the Helmholtz equation), we directly apply li minimization techniques to solve 
the underlying PDEs. This approach is more expensive as it involves a large snapshot space; 
however, in this example, we can not identify a local minimization principle, such as local 
generalized SVD. 

All our numerical results assume the sparsity and we discuss this assumption for the snapshot 
spaces. Moreover, we discuss the computational savings provided by our approach. The sparse 
solution allows a fast evaluation of stiffness matrices and downscaling the solution to the fine 
grid since the reduced dimensional solution representation is sparse in terms of local snapshot 
vectors. Numerical results are presented, which show the convergence of the proposed method 
and the sparsity of the solution. 
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1 Introduction 


1.1 Multiscale problems and the problem of sparsity 

Simulations of multiscale problems are expensive and, typically, require some type of a model 
reduction. Our approaches seek adaptive reduced-order models, locally in space, and construct 
multiscale basis functions in each coarse region to represent the solution space. These approaches 
share common concepts with homogenization and upscaling methods liaElEailSlEHlEolESlEll, 
where local effective properties are constructed. In contrast, in multiscale methods, local multiscale 
basis functions IMlEHlEaEaiinillSlIISlEaillllZlEllllHlEl] are constructed to represent 
the solution space. These basis functions are typically constructed in the snapshot spaces [20] . 
In this paper, we investigate cases, when the basis functions are sparse in the snapshot space. 
We discuss several examples, which include multiscale parameter-dependent problems and the 
Helmholtz equation. The parameter-dependent multiscale problems are motivated by stochastic 
problems, where the parameter is used to describe the uncertainties. 

1.2 Sparse GMsFEM Concepts 

In this paper, we use the GMsFEM framework f [?n[^[^l^ l^ [T^[T^[TTl[T0l[^ ) and investigate 
the sparsity within GMsFEM snapshots. To illustrate the main idea of our approach, we consider 

Lu = f, 

where L is a differential operator. For example, in the paper, we consider parameter-dependent 
heterogeneous ffows, Lu — —div{K,{x]fi)'Vu)^ and the Helmholtz equation, Lu — —div{tv{x)'Vu) — 
Vt^n{x)u. The main idea of GMsFEM is to construct a snapshot space and identify a subspace, 
called the offline or online space depending whether the problem is parameter-dependent. This 
subspace is used to solve the underlying problem at a reduced cost. The snapshot and online 
spaces are constructed in each coarse element (see next section for more precise definitions), where 
a coarse element is a region, which is much larger than the characteristic fine-length scale (see 
Figure]^. For each coarse region, Tj, we construct snapshot vectors, {'0^} (here i is the numbering 
of the snapshot functions), that represent the local solution space. We denote the snapshot space 

Vsn^p = SpanJ?/^^}, Esnap = 

In GMsFEM, the online spaces are constructed using the elements of local snapshot functions. In 
many examples, the snapshot space can be large and the online space can be a sparse subspace of 
the snapshot space. The objective of this paper is to investigate these cases. 

1.3 Snapshot spaces 

The snapshot spaces play an important role in the GMsFEM. They are designed to capture the 
solution space locally and are used to preserve some features of the solution space, e.g., mass 
conservation. Typical snapshot spaces consist of local solutions constructed using some sets of 
boundary conditions or right hand sides. With an appropriate choice of snapshot spaces (e.g., 
using oversampling EH), one can improve the convergence of GMsFEM substantially. 

To convey the concept of snapshot spaces, we present some examples. We will consider two 
examples discussed above. We start with a simplified example related to the parameter-dependent 


2 


case, Lqu — —div{K,{x] fi — 0)Vi/), i.e., the problem without a parameter. In each coarse-grid block 
Tj (see the left plot in Figure]^, we consider a local solution 

Lo(V’0 = 0 in Tj (1) 

subject to some boundary conditions, where these boundary conditions play an important role in 
defining snapshot functions. One option is to choose all possible boundary conditions considering 
all unit vectors on the boundary of Tj. More precisely, {x) = 5i{x) on dxj^ where 5i{x) is 1 at the 
node i and zero elsewhere. The computations of these snapshot functions are expensive. Instead, we 
use the boundary conditions, which are randomly distributed numbers on the fine-grid nodes of the 
boundary dxj (see the left plot in Figure]^. The random boundary conditions allow extracting the 
essential information provided we choose several more snapshot vectors than the number of modes, 
we would like to use. For parameter-dependent problems^ the snapshot vectors are defined as above 
Q for each pre-selected value of /i^. For example, for one-dimensional ease rj = [xj,Xj+i], for 
parameter-independent problem, the snapshot space in tj consists of two solutions and 
such that ^'^{xi) = 6ni^ n — j + I — j, j -h 1, where 6ji is the Kronecker symbol and 
(^ = + 1) is a solution of ^ (a^(x;/ i = 0 )^ 7 /^^) = 0 in tj (see Figure for illustration). For 

parameter-dependent problems^ the snapshot vectors in rj — [xj^Xj^i] are the solutions of 

^( Kte , O ^*)=0 1nr,, 

= ^nU = J 5 J + I 5 ^ + 1- For multi-dimensional examples, we can construct the 

snapshots similarly for each Tj and for using different boundary conditions and use one index to 
represent the snapshot vectors as For the second example, Lu — —div{n{x)Vu) — Q‘^n{x)u^ we 
choose the snapshot vectors to be functions for a set of pre-defined values of ki on a unit 

circle (see the right plot in Figure]^. 



Figure 1: Illustration of snapshot concepts in one dimensional example. 
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The multiscale basis functions are constructed in the snapshot space. Our earlier approaches 
seek a small dimensional subspace of the snapshot space by performing a local spectral decom¬ 
position (based on analysis). However, these approaches use all snapshot vectors when seeking 
multiscale basis functions. In a number of applications, the solution is sparse in the snapshot 
space. I.e., in the expansion 

hj 

many coefficients Cij are zeros. In this case, one can save computational effort by employing sparsity 
techniques. In this paper, our main goal is to discuss how GMsFEM can be designed if the solution 
is sparse in the snapshot space. We describe two classes of approaches and present a framework for 
constructing sparse GMsFEM. 

The main challenge in these applications is to construct a snapshot space, where the solution is 
sparse. In our first example, this can be achieved, because an online parameter value fi can be close 
to some of the pre-selected offline values of /i’s, and thus, the multiscale basis functions (and the 
solution) can have a sparse representation in the snapshot space. In our second example, we select 
cases where the solution u contains only a few snapshot vectors corresponding to directions ki. We 
note that if the snapshot space is not chosen carefully, one may not have the sparsity. In general, 
there can be many other examples and our goal is to show how local multiscale model reduction 
techniques can be used for such problems. 



Figure 2: Illustration of snapshot concepts. Left: For the first problem; Right: For the second 
problem. 


1.4 Approaches for identifying sparse solutions in snapshot spaces 

The main goal of this paper is to discuss how to explore sparsity ideas within GMsFEM for con¬ 
structing local multiscale basis functions. We consider two distinct cases. 

• First approach: “Local-Sparse Snapshot Subspace Approach”. Determining the online sparse 
space locally via local spectral sparse decomposition in the snapshot space (motivated by 
parameter-dependent problems). 

• Second approach: “Sparse Snapshot Subspace Approach”. Determining the online space 
globally via a global solve (motivated by using plane wave snapshot vectors and the Helmholtz 
equation). 
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See Figure for illustration. We use sparsity techniques (e.g., [SI El 0 ESI HU) to identify local 
multiscale basis functions and solve the global problem. 

In above approaches, the snapshot functions can be linearly dependent. In fact, in general, we 
would like to have a large snapshot space that can contain a sparse representation of the solution. 
In both approaches formulated above, the linear dependency is removed. In the first approach, 
it is removed by sparse POD. We note that in the original GMsFEM approach [2Q], POD across 
all snapshot functions is used to remove linear dependence. However, this can result in a loss of 
sparsity, i.e., the solution may contain many nonzero coefficients when represented in the snapshot 
space. Thus, for sparsity, it is important to avoid a POD step across all snapshot vectors. In our 
first example, we need to avoid using POD for all /i’s. For this reason, we design a special sparse 
POD method using randomized snapshot functions. It both eliminates linearly dependent snapshot 
functions and identifies a sparse solution space. In the second method, li minimization can be used 
for all snapshots, even if they are linearly dependent. By adding more snapshot vectors, we hope 
to identify a sparse representation of the solution. The second method will eliminate the linear 
dependence and identify a sparse solution. Note that in our example, the snapshot vectors are 
linearly independent. 



Figure 3: Illustration of our approaches. 

For the first case, we consider parameter-dependent elliptic equations of the form 

— div(/^(x; jj) Vu) = / inD, (2) 

where ^ on dD^ and /i is a parameter. Some of the existing approaches for parameter- 

dependent problems closely related to the proposed approaches include reduced basis techniques 
mu- In these approaches, the reduced order model is constructed via a greedy algorithm. In 
the proposed approach, we attempt to approximate the solution space locally in each coarse block 
using li minimization. Local multiscale basis functions are constructed using GMsFEM. In previous 
approaches, we attempted to compress the local solutions corresponding to some pre-selected values, 
/i, in the offline stage. This can lead to large dimensional offline spaces. In this paper, we propose 
an approach to compute eigenvectors using li minimization based on randomized snapshots and 
oversampling. This provides an efficient approach to identify sparse eigenvector representation in 
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GMsFEM. The proposed approach gives a sparse representation of the multiscale basis functions 
in terms of the snapshot space vectors and has several advantages. (1) It allows quickly assembling 
of the stiffness matrix in the online space since it involves a few elements of the snapshot space. 
(2) We can downscale the solution to the fine grid much faster using sparse representation. (3) 
It avoids a POD based step proposed in the original GMsFEM formulation (see m), which is 
performed across all /i’s and which can result to a large dimensional representation of the online 
multiscale basis functions in terms of snapshot functions. 

In the second example, we consider the Helmholtz equation 

— div(A^(x) Vu) — Q‘^n{x)u = / inD, (3) 

where D is the frequency. We will use plane waves as the snapshot vectors. In the computational 
examples considered in this paper, the solution has a few dominant propagating directions, and 
is, therefore, spanned by only a few plane waves. This observation leads to the solution sparsity 
in our snapshot space. However, the choice of local spectral decomposition is not available for 
determining these dominant directions and we study using li minimization directly in the space of 
snapshot vectors. We consider not very large snapshot spaces, and snapshot vectors having closed 
form formulas for constant media properties (/^(x) and n{x)). For the Helmholtz equation 
with low frequencies, we can expect a sparsity in our examples, which we exploit. Thus, sparsity 
techniques are the natural methodologies in these situations. In general, we can also consider the 
frequency to be a parameter, and apply similar techniques used for the first example. 

1.5 Summary of numerical results 

Two test cases for the first approach are presented, where we consider parametrized conductivity 
fields. In the first case, the conductivity is parametrized as an affine combination of two heteroge¬ 
neous conductivity fields. In the second case, we use a nonlinear parameter dependence. In partic¬ 
ular, we consider an initial conductivity field with a channel and inclusions. The parametrization is 
introduced such that these high-conductivity features spatially move within the domain. This is a 
more challenging example because high-conductivity features appear in many parts of the domain. 
Numerical results show that our approach provides an accurate approximation of the solution using 
a few degrees of freedom and the solution is sparse in an appropriate snapshot space. 

Numerical results for the second approach involve solving the Helmholtz equation in media with 
two isolated heterogeneous inclusions. We consider a domain with two distinct properties, where 
plane wave solutions can provide a good approximation. In this case, the solution is spanned by 
only a few plane waves, and is, therefore, sparse in the space of plane waves with many propagating 
directions. However, in general, we do not know which plane wave directions are dominant and our 
algorithm identifies these directions by using li minimization. 

The paper is organized in the following way. In the next section, we present preliminaries and 
discuss GMsFEM, coarse and fine grid concepts. In Section 3, we propose our new construction for 
the online space. Section 4 is devoted to numerical results. In Section 5, we present conclusions. 

2 Preliminaries 

To discretize ([^ or (|^, we let be a usual conforming partition of the computational domain 
D into finite elements (triangles, quadrilaterals, tetrahedrals, etc.) and denotes all the edges 
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in the coarse mesh . We refer to this partition as the coarse grid and assume that each coarse 
subregion is partitioned into a connected union of fine-grid blocks. The fine-grid partition will be 
denoted by We use {xi}^^ (where Ny denotes the number of coarse nodes) to denote the 
vertices of the coarse mesh and define the neighborhood of the node Xi by 

oJi = \J{Kj G r^; Xi e Kj}. (4) 

See Figure for an illustration of neighborhoods and elements subordinated to the coarse dis¬ 
cretization. We emphasize the use of oji to denote a coarse neighborhood, and K to denote a coarse 
element throughout the paper. 
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Figure 4: Illustration of a coarse neighborhood and coarse element 

Next, we briefly outline the global coupling and the role of coarse basis functions for the re¬ 
spective formulations that we consider. For the discontinuous Galerkin (DG) formulation, we use 
a coarse element K as the support for basis functions, and for the continuous Galerkin (CG) for¬ 
mulation, we use (jJi as the support of basis functions. In turn, throughout this chapter, we use the 
notation 

/ Ui for CG . . 

for DG 

when referring to a coarse region where respective local computations are performed (see Figure]^. 
To further motivate the coarse basis construction, we offer a brief outline of the global coupling. In 
particular, we note that our approach will employ multiple basis functions per coarse neighborhood. 
Both CG and DG solutions will be sought as fi) — ^ /i), where ji) are 

the basis functions (without loss of generality, we write basis functions as parameter-dependent). 
Once the basis functions are identified, the global coupling is given through the variational form 

a*) = (/> v), for all v G (6) 

where is used to denote the space formed by those basis functions and a£)G/CG i® a bilinear 

form which will be defined later on. Throughout, for the convenience, we use the same notations 
for the discrete and continuous representations of spatial fields. 
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3 Sparse GMsFEM 


In this section, we will give the detailed constructions of our sparse GMsFEM. We start with an 
outline of the approach. 

3.1 Outline 

In this section, we present an outline of the algorithm. Assume that the snapshot space is l^nap ~ 
Span{'0^^^^} for a generic element r. We assume that the solution is sparse in this snapshot space 
and consider two approaches (see Figurefor illustration). Throughout the paper, we will assume 
that the solution is sparse in the snapshot space. We will discuss the assumption on the sparsity 
later in Section [5l 

General outline of the sparse GMsFEM: 

1 . Coarse grid generation. 

2 . Construction of snapshot space, where the solution is sparse. 

3. — First approach. Local-Sparse Snapshot Subspace Approach. Seek a subspace of the snap¬ 

shot space and construct multiscale basis functions that are sparse in the snapshot space. 

— Second approach. Sparse Snapshot Space Approach. Solve for the sparse solution in the 
snapshot space directly within a global formulation. 

In the first approach, we perform local calculations to identify multiscale basis functions that 
are sparse in the snapshot space. Here, we will use approaches similar to sparse POD. Then, the 
global problem is solved in the space of multiscale basis functions. The resulting solution is sparse. 

In the second approach, we apply directly sparse solution techniques and find the solution that 
is sparse in the snapshot space. This approach is more expensive since it uses a large snapshot 
space. However, in some examples, we can not identify local basis functions in the offline stage and 
such approaches give sparse solutions in the online stage. 

3.2 First approach. Local-Sparse Snapshot Subspace Approach 

We first give a general idea of this approach. We consider a local snapshot space t^nap ~ 
Span{'0^^^^}. In the local snapshot space, we seek multiscale basis functions that are sparse 

in the local snapshot space and which have smallest energies (similar to sparse POD). For example, 
following to |13] , we can consider 

min — ||T|L + Tr(T^A 5 (/i)T), s.t. = I, 

where n is the dimension of the fine-scale space, is the number of online basis functions and 
As{/jl) is the stiffness matrix formed in the snapshot space. Our proposed local approach avoids 
expensive direct local eigenvalue calculations and uses randomized snapshots. This approach will 
be used for problems where one knows a local basis construction principle. The latter typically 
uses some local generalized eigenvalue problems in the snapshot spaces mi. Once multiscale basis 
functions (that are sparse in the snapshot space) are constructed, we solve the problem on a coarse 
grid. 


Here, we will consider the parameter-dependent problem (§. We can consider the computation 
of the parameter-dependent coarse space as an online procedure. In the latter, our objective will 
be solving many problems for a given value of the parameter with different boundary conditions 
and the right hand sides. 


3.2.1 Snapshot space 

We first construct a snapshot space (for a generic r). Construction of the snapshot space 

involves solving the local problems for various choices of input parameters, and we describe it 
below. We generate snapshots using random boundary conditions by solving a small number of 
local problems imposed with random boundary conditions. 


—div(/^(x; = 0 in 

/ r.rsnap -u 

^ = r/ on dr^. 


(7) 


where ri are independent identically distributed (i.i.d.) standard Gaussian random vectors on the 
fine-grid nodes of the boundary, / = !,••• , L. /ij (j = 1,..., J) is a specified set of fixed parameter 
values, and J denotes the number of parameter values we choose. Here, is an oversampled 
region shown in the Figure as uf or for conforming Galerkin formulation or discontinuous 
Galerkin formulation. 

The space generated by is a subspace of the space generated by all local snapshots 

A; = !,••• ,A^, and N denotes the number of boundary nodes. Denote = 

brsnap^ • • • , and "P, • • • ,^”]. Therefore, for each parameter /i„ there 

exists a randomized matrix TZ with rows composed by the random boundary vectors ri (as shown 
in Equation Q), such that. 


^r,rsnap _ .^^r,snap 


( 8 ) 


Now, we are ready to present the local snapshot space as follows. 


VLp = : 1 < j < J and 1 < I < L}, 

for each coarse subdomain r. 


Remark 3.1. Note that we impose the same random veetors for the loeal snapshot ealeulation in 
Equation 0 for /ij, j = 1, • • • , J in order to obtain a sparse online spaee. 


3.2.2 Sparse local space calculations 

For a given input parameter /i, we next construct the associated coarse space for eaeh pi value 

on each coarse subdomain r. In principle, we want this to be a small dimensional subspace of the 
snapshot space for computational efficiency. The coarse space will be used within the finite element 
framework to solve the original global problem, where a continuous or discontinuous Galerkin 
coupling of the multiscale basis functions is used to compute the global solution. In particular, 
we seek a subspace of the snapshot space f^nap such that it can approximate any element of 
the snapshot space in an appropriate sense. For the convenience of the presentation, we denote 
^snap ~ • • • , Similar as the generation of local snapshot space, we obtain the local 

sparse online basis via a local problem solved by li optimization in the space of the corresponding 
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local snapshot space. Here, we will use a smaller subspace of V^nap space constructed 

through multiplication of a random matrix r^andom* Denote T^andom ^ matrix of size (L x J)by q 
{q << {Lx J)) with rows of i.i.d standard Gaussian random vectors. Then the test space is defined 

oq \[/T rpr 

^ snap random ‘ 

Specifically, the local problem is arranged as follows. Find such that, = ^Lap^^ 

Ui — argmin- \\U\\-^ subject to Ac{jJ^)Ui = F/. (9) 

Here, Ad^i) = {'^LapT^ndomf and Fi = Ri with being the local 

stiffness matrix and Ri being the right hand side for the local problem with Dirichlet boundary 
condition ri. Namely, we are solving the following local problems in l^^ap with the li minimized 
coefficient vector Ui for the testing space equal to ^Lap^random’ 


—div(/^(x; /i) 


,on\ 


= 0 in 


= ri on dr^. 


( 10 ) 


Later, we will briefly introduce the algorithm to solve Equation Q. Note that we impose the same 
random vectors as the boundary conditions in Equation Q to guarantee the sparse solution in 
Equation (10). 


We can obtain the local online snapshot functions on the target domain r by restricting the 


solution of the local problem, to r (which is denoted by 

the local online snapshot space as follows. 


Now we are ready to present 


Vn = Span{V’[’°" : 1 < ^ < i^}, 

for each coarse subdomain r. Then we denote • 5 

Next we will select the dominant modes from by the following eigenvalue problem, 

( 11 ) 


where 

5---(m) = [s^dl^Un] = 


and i^{x] fi) is now parameter dependent. To generate the coarse space we then choose the smallest 

eigenvalues from Equation 0 and form the corresponding eigenvectors in by setting 
L 

^ A: = 1,..., where are the coordinates of the vector 

j=i 


Above, we presented an algorithm for solving an eigenvalue problem in the space which has a 
sparse representation in the snapshot space. One can attempt to obtain sparse spectral basis from 

([ansT]). Here, we 


Equation (10) using /i-minimization method in the local snapshot space 
can use the algorithm proposed in m- The general basis pursuit problem is as follows 


rr 
snap 


min llxIL + u \\Cx — f{x)\\ , 
xeR^ 


( 12 ) 
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where / G C G and m « n. We refer to [U] for the Bregman algorithm to solve 


(12). Instead of solving Equation (10) in the online stage, we can generate the sparse spectral basis 


directly following the algorithm in m- The first sparse spectral basis can be solved by 


min 1 ll^lli + s.t. = I. 

1/ 


(13) 


3.2.3 Global coupling 

The multiscale basis functions constructed above can be coupled via DG or CG formulation. Below, 
we present DG approach (similar results are observed when CG approach is used). One uses the 
discontinuous Galerkin (DG) approach (see also [l3l[T]) to couple multiscale basis functions. This 
may avoid the use of the partition of unity functions; however, a global formulation needs to be 
chosen carefully. The global formulation is given by 


aDG{uW^,v) = {f,v), yvEVon, 


(14) 


where the bilinear form is defined 


as 


aDG{u,v) ^ aHiu,v) - ^ J ({kVw • n£}[n] + {kVv • n^jM) + E ^ / '^MH (15) 


Ees^ 


Ees^ 


with 


aniu^v) = <^^('^ 5 '^)? a^{u^v) = / kVu’Vv^ 


(16) 


K&Th 


where 7 > 0 is a penalty parameter, is a fixed unit normal vector defined on the coarse edge 

E G . Note that, in rtl5L the average and the jump operators are defined in the classical way. 


Specifically, consider an interior coarse edge E G and let and K be the two coarse grid 
blocks sharing the edge E. For a piecewise smooth function G, we define 


{g}E-{g+ + g- 


|G] = G+-G-, on E, 


where = G|x+ and G~ = and we assume that the normal vector is pointing from 

to K~. Moreover, on the edge E, we define K = {i^k+ + ^k-)/ 2 where is the maximum value 
of K. over K^. For a coarse edge E lying on the boundary dD^ we define 

{G} = |G] = G, and 'R — k^k on F, 


where we always assume that tie is pointing outside of D. We note that the DG coupling (14) is 
the classical interior penalty discontinuous Galerkin (IPDG) method [l3] with our multiscale basis 
functions as the approximation space. 

We can obtain the discontinuous Galerkin spectral multiscale space as 


Vn'^(M) = Span{^“ 


K,on 






(17) 


We can obtain an operator matrix constructed by the basis functions of V^n^(/i). We denote 
the matrix as $0 where $0 = [ 0 ?^? 


, 1 /)^^]. Recall that Nc denotes the total number of coarse 
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basis functions. Solving the problem (§ in the coarse space V^n^(/i) using the DG formulation 
described in Equation (14) is equivalent to seeking /i) = fi) G such that 


= {f,v) for all V G ( 18 ) 

where (jP^{u^v] jh) and f{v) are defined in Equation ( [T^ . We can obtain a coarse system 

AoU^^ = Fo, (19) 

where denotes the discrete coarse DG solution, and 

Ao{ii) = i?o ^(/i)i?o and Fq = RqF^ 

where A{fi) and F are the standard, fine-scale stiffness matrix and forcing vector corresponding to 
the form in Equation (18). After solving the coarse system, we can use the operator matrix Rq to 
obtain the fine-scale solution in the form of RqU^^. 


3.2.4 Computational cost 

In this section, we discuss the computational cost. For this, we assume that we have chosen J 
parameters, and L boundary conditions ri,...,rL, for constructing the snapshot space. 

Then, the cost for snapshot calculations will be the same as solving J x L local problems for 
randomized snapshots. Next, we compute the cost of solving Lon online randomized snapshots. 
Each online snapshot calculation requires solving li minimization with a constraint involving q x 
(L X J) matrix (see ©)• The cost of eigenvalue computation with Lon snapshots is considered to 
be small as it involves a small eigenvalue problem of the size Lon x Lon- The online cost is mainly 
due to solving (1) solving Lon online randomized snapshots (2) solving a global problem on a coarse 
grid. The cost of solving a global problem is small if the solution has a sparse representation. As 
for the cost of solving online randomized snapshots, this will be small compared to solving local 
problems in the online stage if local problems have a high resolution. Moreover, as we pointed out 
earlier that the proposed approach allows a fast assembly of the stiffness matrix in the online space 
since it involves a few elements of the snapshot space. It also avoids using all snapshot vectors, 
which can result to a large dimensional representation of the online multiscale basis functions. 


3.3 Second approach. Sparse Snapshot Subspace Approach 

In this approach, we will use an appropriate snapshot space to compute the sparse solution directly. 
Again, we consider a local snapshot space V^nap ~ SpanlT/^J’^^^^}. In some applications, we may 
not be able to reduce the dimension of the multiscale space locally. In this case, we can use sparsity 
techniques directly in the global snapshot space to compute the solution. The procedure can be 
more expensive; however, can yield more accurate solutions. More precisely, we seek the solution 
in the global snapshot space of l^nap = Spanl?/^^^^^} using li minimization with testing space, 
Vtest^ spanned by the random combination of snapshot basis functions. This is equivalent to find 
Urns = Y.i e Vnap where 

= argmin ||;7||;^ subject to aDG(E] '^) = (/’ ^ ^test- (20) 
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In the following section, we will take the Helmholtz problem as an example and discuss the 
procedure of using this approach to compute a sparse solution. More precisely, we consider the 
following problem: find u such that 

—V • — VL^n{x)u = /, in D 

with the Dirichlet boundary condition u\qd = g. 

3.3.1 Snapshot and test space 

In this section, we present the construction of the snapshot space l^nap and the test space Vtest- 
Since we are solving the Helmholtz equation with a fixed frequency f], we can assume the solution, 
can be written as a linear combination of plane waves, namely, u = • Therefore, we 

consider our snapshot basis to be some plane waves in each coarse block K G that is, 

V^nap = Spnn{^rnj • ^ < m < Nd and 1 <j < M}, 


where 


(^) 


0 otherwise 


( 21 ) 


with kjn = (sin(7rm/A^^),cos(7rm/A^^)) , Nd as the number of propagating directions and M the 
number of coarse blocks. We note that the plane wave basis is not new in solving the Helmholtz 
equation and it was used in a number of papers (see [aaiiHiiiHiiiaES] and the references therein). 

Next, we will show the construction of the test space. We consider as a collection of 

i.i.d. standard Gaussian random vectors with Nt « Nd- Then the testing space l^est is defined by 


Nd 

^est — — E 1 < I < Nt and 1 < j < M }. 

m=l 


Since Nt « Nd^ we have a test space with dimension much smaller the snapshot space (dim(Vi;est) = 
NtM « NdM = dim(Hsnap)). 


3.3.2 Sparse solution in the snapshot space 

After constructing the snapshot and test space, we can couple the system globally by IPDG method. 
That is, we find u^s ^ ^nap such that 

<^DG('^ms5'^) — (Z?'^)? ^ ^est? (2^) 


where 


aBG{u,v) 


aH{u,v) 


“ / ({i^Nu-nE}lvj + {KVv-nE}luj) + ^MN 


(23) 


with 


aniu^v) 


a^{u,v), 

kgTh 



kVu ■ Vv — 



(24) 
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where 7 > 0 is a penalty parameter, is a fixed unit normal vector defined on the coarse edge 
E G . Moreover, {•} and [•] are the average and jump operators defined before. 

As the dimension of test space is smaller that the dimension of snapshot space, the Equation 


(|22|) does not have a unique solution. To seek for a sparse solution, we will solve a li minimization 

, more precisely, we will find u 


problem subject to Equation 


- '^i,j 


= argmin{||{7||/i} subject to aDG{Uu4>l,j,v) = (/,u), Vu € Ftest 


G f^nap such that 


(25) 


3.3.3 Cost of computations 

In this section, we discuss the computational cost associated with the second approach. In this 
example, the cost for snapshot (plane waves) calculations is cheap as they are analytically de¬ 
scribed. In addition, the linear system in ([2^ has dimension NfM x which is a highly 


under-determined system as Ni « Thus, one can solve the li minimization problem (25) effi¬ 


ciently by using, for example, the Bregman’s method m- Moreover, if an adaptivity can be used 
and the problem requires very few snapshot vectors or very few test vectors in many regions excepts 
a few coarse regions, this will increase the efficiency of the proposed approach. In this case, if we 


denote by and the number of test and snapshot vectors in the region i (using adaptivity). 


ri'^) 


then the linear system in (25) has dimension 


is not large, one can gain computational efficiency. We note that this approach is more 
expensive compared to the first approach, where we perform the sparsity calculations at the local 
coarse-grid level. 


/"(d ‘c 


Ar(i) 


Ar(b 


Thus, if ZZi 




or 


4 Numerical results 

4.1 First Approach. Local-Sparse Snapshot Subspace Approach 

In this section, we will present some numerical examples by using the first approach to compare the 
sparse multiscale solution. We consider the domain D = [0,1]^. The coarse mesh size H is 1/10 
and each coarse grid block is subdivided into a 10 x 10 grid, therefore, the fine mesh size h — 1 / 100 ; 

Example 1 

Setup. In our first example, we will consider the source function f — 1 and the medium parameter 
A^(/i) = (1 — /i)/^i + where ni and are shown in FigureWe choose the offline values of /x, 
jdi = 0.2,0.4, 0.6, 0.8, for computing the online snapshot space as discussed above. 

Discussions of numerical results. In Table[^ we show the convergence history of our method 
for ji = 0.5, where we defineThe fine-grid solution and the numerical 
solution are shown in Figure]^ First, we note that there is an irreducible error due to the use 
of the snapshot space, which consists of harmonic functions. This error is of order of the coarse 
mesh size and this is the reason, the error decay is slow (below 10 %) as we increase the dimension. 
In these problems, because our selected fi = 0.5 is near to /x = 0.4 and /x = 0.6, we observe that 
the sparsity is close to 50 %, i.e., we only use snapshots corresponding to nearby values of /x. We 
observe that when the snapshot space dimension is 9600 (i.e., 24 randomized solutions per coarse 
block and per each value of /x^), the nonzero coefficients in the expansion of basis functions (over 
the whole domain) in terms of 9600 snapshot vectors are 4850. The optimal expansion for this case 
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Figure 5: Decomposition of permeability field 


corresponds when ji — 0.5 is selected for snapshot construction and in this case, the number of 
nonzero coefficients is 2400 (24 per coarse region). We note that if we consider small dimensional 
online spaces and the full snapshot space, then the sparsity is very small. For example, if we use 
12 randomized solutions per coarse block and per each value of jii for identifying multiscale basis 
functions per each coarse region, then, we will be using only 1/2 of the snapshot vectors and thus, 
the sparsity (the number of nonzero coefficients of the solution in the snapshot space) will be 25 
%. As we observe that our numerical examples identify appropriate sparsity of the solution space. 
We expect a more significant gain in the sparsity if more parameter values are used. 

Why to expect a sparsity. Next, we briefly describe why to expect a sparsity in this problem. 
Because the snapshot space consists of local problems corresponding to multiple values of /i, we 
expect that for an online value of /i, we will have a local (in K) coefficient /^(x; /i), which is similar 
to one of the snapshot solutions. Thus, it is more advantageous to use li minimization techniques, 
which will select a small dimensional subspace of the snapshot space corresponding to the coefficient 
that is close to a^(x; ji) with the online value of /i. Such situations may occur in various applications. 
Moreover, in these examples, it is more advantageous to use local spectral decomposition and avoid 
a large-scale h minimization problem. 


Table 1: Convergence history of the DGMsFEM using oversampling Harmonic basis. The fine-scale 
dimension is 12100. The full snapshot space dimension is 9600. 


dim(yon) 

1 '^msl (^) 

L'\D) 

H^iD) 

400 

15.05 

31.84 

600 

2.89 

13.71 

800 

1.22 

10.20 

1000 

1.12 

9.83 

dim( V^nap) 

1.07 

9.59 
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Figure 6: Left: Fine grid solution Right: Numerical solution (8 basis). 


Example 2 

Setup. In our second example, we consider the source function / to be the same as in the previous 
example. As for the medium parameter, we use a nonlinear function where the high permeability 
channel and inclusions move as we change the parameter. The expression for the medium parameter 
is 

K(/i) = Ki(a; + /i,y) 

In Figure]^ we show four different values of n, {fj, = 0,0.15,0.3,0.45), that is used to construct the 
snapshot space. This is a complicated case as the high-conductivity region is not fixed. 

Discussions of numerical results. In Table we show the convergence history of our 
method. The fine grid solution and the numerical solution are shown in Figure with fi = 0.14. As 
we see from this table that the error is larger compared to the previous case. The error decreases as 
we increase the dimension of the space. However, due to the fact that we do not span all (or many) 
parameter values, the decay is slow. As for the sparsity, we achieve a better sparsity compared to 
the previous example because the online value of /i = 0.14 is close to one of selected offline values 
fi = 0.15. In fact, we observe that when the snapshot space dimension is 9600, as before, (i.e., 
24 randomized solutions per coarse block and per each value of /i^), the nonzero coefficients in the 
expansion of basis functions (over the whole domain) in terms of 9600 snapshot vectors are 3700. 
The optimal expansion for this case (i.e., the case with 24 randomized solutions per coarse block) 
corresponds when /i = 0.15 is selected for snapshot construction and in this case, the number of 
nonzero coefficients is 2400 (24 per coarse region). If we consider small dimensional online spaces 
and the full snapshot space, then the sparsity is very small, as before. For example, if we use 
6 randomized solutions per coarse block and per each value of fii for identifying multiscale basis 
functions per each coarse region, then, we will be using only 1/4 of the snapshot vectors and thus, 
the sparsity will be 9.5 %. As we observe that our numerical examples identify appropriate sparsity 
of the solution space. Again, we expect a more significant gain in the sparsity if more parameter 
values are used. 

Remark 4.1. We have implemented CG-GMsFEM using speetral basis approaeh and observed 
similar results. 
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dim(yon) 

ll«-«°"ll (%) 

L^D) 

mD) 

800 

13.50 

31.06 

1000 

12.02 

29.45 

1200 

9.72 

26.66 

1400 

7.97 

24.13 

dim(14nap) 

5.78 

20.27 


Table 2: Convergence history of the DGMsFEM using oversampling Harmonic basis. The fine-scale 
dimension is 12100. The full snapshot space dimension is 9600. 
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Figure 7: medium parameter. Top-Left: Top-Right: Bottom-Left: Bottom- 

Right: 

4.2 Second Approach. Sparse Snapshot Subspace Approach 

In this section, we will show the numerical example by using second approach to directly calculate 
the sparse multiscale solution by li minimization. 

Example 1 

Setup. In this example, we consider the domain D = [0,1]^ that is partitioned into the coarse grid 
with grid size iL = 1/8 and each coarse block is subdivided into 16 x 16 fine square blocks with 

length h — H/IQ. Therefore, the fine mesh size h — We consider Vt — 2^ n = 1 and n{x) is 
shown in Figure]^ We consider a zero source function with Dirichlet boundary condition given 
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Figure 8: Left: Fine grid solution Right: Numerical solution (14 basis). 


by ^ = e where k — (sin(—), cos(—)). 

Discussions of numerical results. We will compare our result with the reference solution, 
which is calculated on the fine grid and shown in Figure [T0| Notice that, within each coarse grid 
block, the reference solution has few dominant propagating directions, which suggests sparsity of 
the solution in the snapshot space. In this case, the snapshot space is spanned by local plane waves 
with dimension dimCl^nap) = 1280, as defined in ( [2T| ) with /c^’s distributed uniformly. The snapshot 
solution (i.e., if we use all snapshot vectors) has 1.63% relative error with respect to the fine-scale 
solution. We compare the solutions in Figure [TT| As we observe, the snapshot solution is accurate. 
Next, we calculate the sparse solution by varying the dimension of the test space. The latter defines 
a sparse solution in the subspace of the test space. The numerical solution calculated with 4 test 
basis per coarse grid block is shown in Figure [T^ In Table we show the convergence history 
of the second approach, where \Vu\‘^. As we observe that for low dimensional test 

spaces, the solution is very sparse in the snapshot space (and this sparsity is about the same as the 
test space). We increase the dimension of the test space to achieve a higher accuracy. For example, 
for the solution with the sparsity 408 (i.e., 408 non-zero coefficients in the span of 1280 snapshot 
vectors), we have 1.63 % error. 

Why to expect a sparsity. Next, we briefly describe why to expect a sparsity in this problem. 
The snapshot space consists of local problems corresponding to different directions km (see ([^). In 
this problem, we expect that the solution will consist of plane wave solutions with a few directions. 
By using plane wave solutions, we can identify these few directions. Note that in this example, we 
can not identify local spectral decomposition. 


dim(Viest) X 2 

1 '^msl (9o) 

sparsity of the sol 

L^iD) 

H^{D) 

128 

41.10 

195.38 

128 

256 

21.12 

45.16 

252 

384 

14.62 

32.62 

344 

512 

3.49 

14.40 

408 

dim(14nap) 

1.63 

9.88 

1280 


Table 3: Convergence history of the DGMsFEM to compute sparse multiscale solution directly. 
The fine-scale dimension is 16384. The snapshot space dimension is 1280. 
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Figure 9: Parameter n{x). 
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Figure 10: Reference solution Left: Real part of the solution; Right: Imaginary part of the 
solution. 
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Figure 11: Snapshot solution (i.e., when using all snapshot vectors), Left: Real part; Right: 
Imaginary part. 


Example 2 

Setup. In this example, we consider the domain D = [0,1]^ is partitioned into the coarse grid with 
grid size H — 1/16 and each coarse block is subdivided into 16 x 16 fine square fine block with 

length h — Lf/16, therefore, the fine mesh size h = -. We consider = 8 and n(x) is shown 

in figure The parameter k, source function /, and boundary condition g, are the same as the 
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Figure 12: Numerical solution with 4 test basis functions, Left: Real part of the solution. Right: 
Imaginary part of the solution. 


previous example. Because of higher value of we take the fine grid 2 times finer. 

Discussions of numerical results. We will compare our results with the multiscale approach 
with the reference solution, which calculated on the fine grid and shown in Figure [14| Notice that, 
within each coarse grid block, the reference solution has few dominant propagating directions, which 
suggests sparsity of the solution in the snapshot space. In this case, the snapshot space is spanned 
by local plane waves with dimension dim(4^nap) = 5120, as defined in (21) with /c^’s distributed 
uniformly. The snapshot solution error is 2.44% and it is shown in Figure 15 As we observe the 


snapshot solution is accurate. Next, we calculate the sparse solution by varying the dimension of 
the test space. The latter defines a sparse solution in the subspace of the test space. The numerical 
solution calculated with 4 test basis per coarse grid block is shown in Figure [T^ In Table we 
show the convergence history of the second approach. As we observe that for low dimensional test 
spaces, the solution is very sparse in the snapshot space. We increase the dimension of the test 
space to achieve a higher accuracy. For example, the solution with 1958 nonzero coefficients in the 
snapshot space provides 4.25 % accurate solution in sense. 
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Figure 13: Parameter n{x). 
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Figure 14: Reference solution Left: Real part of the solution; Right: Imaginary part of the 
solution. 



Figure 15: Snapshot solution (i.e., when using all snapshot vectors), Left: Real part; Right: 
Imaginary part. 




Figure 16: Numerical solution with 4 test basis functions, Left: Real part of the solution. Right: 
Imaginary part of the solution. 

5 Conclusions 

5.1 Summary of the results 

In the paper, we develop approaches to identify sparse multiscale basis functions in the snapshot 
space within GMsFEM. The snapshot spaces are constructed in a special way that allows sparsity 
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dim(Viest) X 2 

1 '^ms| (^) 

sparsity of the sol 

L^iD) 

H^{D) 

512 

79.12 

176.91 

500 

1024 

74.07 

103.93 

913 

1536 

40.69 

51.79 

1294 

2048 

17.27 

23.74 

1591 

2560 

4.25 

8.63 

1958 

dim(14nap) 

2.44 

6.18 

5120 


Table 4: Convergence history of the DGMsFEM to compute sparse multiscale solution directly. 
The fine-scale dimension is 66049. The snapshot space dimension is 5120. 


for the solution. We consider two apporaches. In the first approach, local multiscale basis functions 
are constructed, which are sparse in the snapshot space. These multiscale basis functions are con¬ 
structed by identifying dominant modes in the snapshot space using li minimization techniques. 
As for the application, we consider parameter-dependent multiscale problems. In the second ap¬ 
proach, we apply li minimization techniques directly to solve the global problem. This approach is 
more expensive as it directly deals with a large snapshot space. As for the application, we consider 
Helmholtz equations. For both approaches and their respective applications, we present numerical 
results and discuss computational savings. Our numerical examples are simplistic and are designed 
to convey the main idea of the proposed approach. 

5.2 Sparsity assumption 

Both approaches assume that the solution is sparse in the snapshot space. The latter requires 
special snapshot spaces, which can yield this sparsity. For example, for local snapshot vectors 
considered in the paper, this requires identifying boundary conditions for the snapshot solutions, 
which can sparsily represent the solution. This may not be easy in general, though in some examples 
can still be achieved. Besides examples presented in the paper, one can consider scale separation 
cases and use piecewise linear boundary conditions. Whether a general framework for constructing 
these snapshot vectors is possible remains an open question. 

5.3 Adaptivity and online basis functions 

In general, once offline spaces are identified, we can use adaptivity mu and online basis functions 
mm to achieve a small error. The adaptivity is accomplished by identifying the regions with 
large residuals and enriching the spaces in those regions. In our earlier works mu, we have shown 
that one needs to use some special error indicators. For the first approach, we can use the “next” 
eigenvector obtained from local eigenvalue decomposition to construct multiscale basis functions. 
For the second approach, one can increase the test space for additional multiscale basis functions. 

In the regions with largest residuals, we can also use online basis functions to reduce the error 
substantially. Online basis functions are computed locally and identified as the localized basis 
functions which can give a largest reduction in the error. These basis functions involve solving 
local problems with a residual on the right hand side (see m for online basis functions for DG). 

In this paper, we can apply adaptivity and online basis functions as discussed in PEU]. For 
parameter-dependent problems, one can consider identifying online basis functions for a set of /ij’s 
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following the analysis in m- This will give a local eigenvalue problem. Another important problem 
for our future consideration is to identify the values of /xi,..., /xj by adaptivity. 
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